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In this paper we examine the methodology for extracting the Knudsen number (K) and the ratio 
of shear viscosity to entropy density (r//s) developed by Drescher et al [J. The final result for 
r\js turns out to be quite sensitive to Glauber parameters, and particularly the parameter x which 
controls the balance between N par t and N co u. We also explore how alternative formulations of the 
functional relation between the elliptic flow and Knudsen number (K = X/R) impacts the physics 
conclusions, based on Pade approximates. Finally, we extend the calculation to include a limiting 
minimum value on the mean free path proportional to the DeBroglie wavelength. These results 
emphasize the importance of clarifying the initial state used in different calculations, as well as the 
ambiguities inherent in using a transport approach in a strongly-coupled regime. 

PACS numbers: 25.75.Dw 



INTRODUCTION 

The medium created in heavy ion collisions at the 
Relativistic Heavy Ion Collider (RHIC) defied expecta- 
tions by showing a strong collective flow, characteris- 
tic of a perfect fluid 0, 03, El, HI- This behavior is par- 
ticularly striking since even a modest amount of vis- 
cous damping, parametrized as the ratio of shear vis- 
cosity r) to the entropy density s, is predicted to result 
in large deviations from ideal hydrodynamics 0. Sev- 
eral methods, each subject to as-yet uncontrolled sys- 
tematic uncertainties, have been developed to estimate 
the value of (v/s) from the experimental data on collec- 
tive flow [l|, 0] , fluctuations Q , and heavy quark trans- 
port A parallel effort in this direction has been the 
development of hydrodynamic codes in two and three di- 
mensions with a properly causal and stable relativistic 
treatment of viscous effects to second order in the ve- 
locity gradients [IE EI EE EE El ■ while recent work 
has resulted in a greatly improved understanding of the 
various formalisms used by these authors [3 , and direct 
comparison to experimental data has been made for a 
range of values in rj/s, there remains considerable debate 
regarding the sensitivity to initial conditions, equations 
of state, and precise freeze-out conditions. 

These different approaches are in large part motivated 
by the conjecture [161 ] that there is a fundamental bound 
on the viscosity to entropy density ratio: rj/s > l/Air. 
Intriguingly, all methods produce values for rj/s from the 
RHIC fluid that are within factors of 1-4 of the bound. 
This brings into sharp focus the need to understand in 
detail the utility and the limitations of each method. In 



this paper we concentrate on the approach of Drescher 
et al. [1], which provides a useful framework for param- 
eterizing departures from ideal hydrodynamic behavior 
in terms of the Knudsen number K = A/ R, where A is 
the mean free path of particles or quasi-particles in the 
system and R is a characteristic measure of the system 
size. 

It has been realized since the earliest applications of 
hydrodynamics to nuclear collisions (see the first two 
footnotes of Ref. (3) that the determination of K is di- 
rectly related to the viscosity of the fluid. For sufficiently 
small values of K the hydrodynamic limit is reached, but 
larger values imply the presence of a non-trivial mean free 
path, leading to the damping of momentum anisotropies, 
and thus deviations from ideal fluid behavior. In Ref. 1] 
the Knudsen number is directly related to the elliptic flow 
scaled by the initial state eccentricity 1)2 /e, thereby pro- 
viding a convenient link between the experimental data 
and the transport properties of the medium. 



METHODOLOGY 

In this section we review the previous extraction of the 
Knudsen number and the the ratio rj/s as formulated in 
Ref. In doing so, we demonstrate the sensitivity to 
various assumptions at different steps in the calculation. 
Following Bhalerao et al. [3, Ref- Q takes the charac- 
teristic size of the system R as the scale of the strongest 
gradient in the initial matter configuration of the system, 
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estimated by 



R = 



V<2/ 2 



(1) 



This expression for R is used to determine the Knudscn 
number, which in turn is assumed to describe deviations 
of the the elliptic flow parameter (V2) scaled by the ini- 
tial eccentricity e from the corresponding ratio for ideal 
hydrodynamics (v2/e)ih vi a 



?-(?) 



1 



e / ih 1 + K/Kq 



(2) 



where Kq is a parameter estimated to be approximately 
0.7 from varying the particle cross-section in a 3D trans- 
port code [19(, or by comparison to a Monte Carlo solu- 
tion of a two-dimensional Boltzmann equation [20| . The 
parametric dependence on K in Eqn.[2]is not directly mo- 
tivated by an underl ying microscopic theory, but rather 
was proposed in Ref. [l8j to have the correct limits in the 
two extremes of large and small K. In the limit of large 
mean- free path (large K) 1)2 /e ^ V-^> while in the small 
mean- free path (small K) limit v^/t approaches the lim- 
iting ideal fluid or hydrodynamic value {1)2/ f)ih with cor- 
rections linear in K. If the functional form of Eqn. [2] was 
unique (to be discussed below), then it is clear that by 
combining the value of {^)ih calculated using ideal hy- 
drodynamics with an experimental measure of t^/e would 
permit a direct determination of K (albeit with the em- 
bedded uncertainties of hydrodynamic calculations men- 
tioned in the introduction). Even more appealing is the 
possibility developed in Ref. [l[ of directly determining 
both (-)ih and K by fitting the experimental data on 
V2/e as a function of centrality. 

The Knudsen number varies with centrality both di- 
rectly via associated changes in R and indirectly through 
changes in the particle number density n: 
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K 



R 



Rna 



(3) 



where a is the effective inter-particle cross-section. 
Please note that having well-defined values for n, a and 
A implicitly assumes classical ballistic transport. The 
prescription for estimating n from the experimental data 
comes from Ref. [181 



1 dN 



Stt dy 



(4) 



which relies on calculations (for example from Monte 
Carlo Glauber [2l[) to determine the transverse area 
St of the system and the assumption of Bjorken 
expansion[22J] with r as an appropriate proper time in a 
rapidity slice. The characteristic time over which the gra- 
dients develop was assumed in Ref. 18] to be r = R/c s 



where c s is the speed of sound from the equation of state, 
giving: 

n(r = R/c s ) 
and thus (by Eqn. [5}: 

1 



1 c s dN 
St R dy 



(5) 



<7C„ dN 



K St dy 



(G) 



Note that with these assumptions the Knudsen number 
no longer requires knowledge of R. One can then elim- 
inate the Knudsen Number entirely and write Eqn. [5] in 
terms of these new quantities. 
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which indicates that fitting plots of the experimental val- 
ues of V2/e as a function of another 'experimental quan- 
tity T^^p can be used to determine (<rc s ) and {^f)ih- 
This approach is then applied to PHOBOS measurements 
of V2 in Cu+Cu and Au+Au collisions at ^Jsnn = 200 
GeV for unidentified hadrons at pseudo-rapidity rj = 
as a function of the number of participating nucleons 
(Np ar t) [H, HI El- We note here and below that the 
dN / dy that appears in these equations is that for to- 
tal particle (presumably parton) rapidity density, so that 
further (plausible) assumptions are required to express 
this in terms of the experimentally measured charged- 
particle pseudo-rapidity density dN^/dr]. 



DETAILS OF THE INPUT PARAMETERS 

In this section we describe the uncertainties associ- 
ated with the determination of the eccentricity-scaled el- 
liptic flow 1)2/ c and the parton transverse areal density 
-J-4r^- Close examination finds that the PHOBOS data 

St dy 

presented in Fig. 1 of Ref. 1] appears quite different from 
what is nominally the same data presented in Fig. 6 of 



Ref. [24j. This is explained primarily by the fact that 



while the extraction of the Knudsen number is based on 
experimental data, the approach relies on various geo- 
metric quantities (e.g. the area St, the eccentricity e) 
which are not measured directly, but which are only esti- 
mated for each centrality bin. Detailing the differences is 
instructive and points out systematic uncertainties that 
must be addressed before a fully quantitative estimate of 
K can be made. 

One significant source of uncertainty in all such calcu- 
lations is the lack of knowledge about the initial distri- 
bution of energy and matter that is relevant for the cal- 
culation of the initial eccentricity e and the overlap area 
St- In Ref. [H, two initial conditions are presented repre- 
senting different assumptions about the initial state, one 
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from a Monte Carlo Glauber (MCG) calculation and the 
other from a particular Color Glass Condensate (CGC) 
calculation [261 ] . For the MCG case the initial distribu- 
tion of energy is defined by a combination of spatial co- 
ordinates of participating nucleons and binary collisions. 
The relative weights of the contributions follow the "two- 
component" model of Kharzeev and Nardi 27( : 



dN, 



ch 



drj 



1 j>{> 



(1-x) 



N, 



part 



• xN, 



coll 



(8) 



with x = 0.20, the so-called "80:20" mixture. This 
nomenclature is potentially misleading: In central 
Au+Au collisions N co a ~ 6N part , which implies that 
the N part term contributes only 40% of the matter while 
N co u contributes 60%. It is also notable that the cho- 
sen value of x = 0.20 is not commonly used in the lit- 
erature (e.g. Ref. [13]) and the experimental data at 
mid-rapidity at all beam energies is best described by 
x = 0.13 ± OM(stat) ± 0.05(sys) [H. 

Fig. Q] shows a new calculation of the eccentricity values 
using a modified version of the PHOBOS MC code [!§] 
where we have incorporated the two-component model. 
Note that we use the standard nuclear parameters and no 
nucleon-nucleon hard core potential (d = 0). We calcu- 
late the eccentricity values from the Monte Carlo event- 
by-event by the following equation: 



tpart 



(9) 



The notation e par t will be used even for x ^ 0, i.e. 
when the weighting is not just for the participants. For 
comparisons with PHOBOS data, the second cumulant 

£part{2} = ((part) i s used, since the event plane method 
is effectively a two particle correlation [iEj ]. 

Fig.fflshows the value of e par t{2} for Au+Au collisions 
as a function of N part for four values of x = 0.0, 0.13, 

0. 20, and 1.00. The lower panel of Fig.[T]shows the ratio 

of each case relative to the x = 0.00 case. Results from 
I — I 

Ref. [1( for x = 0.20 are also shown. We find that our 
results for x — 0.20 are systematically larger than their 
results. This discrepancy has been traced to the fact 
that the analysis in Ref. [l| uses a non-standard form for 
Eqn. with a weight of one (as opposed to 1/2) for the 
N P art term. This difference then gives a larger weighting 
for the spatial distribution of participants (i.e. lower x 
value), though there is no mapping onto a different exact 
x value given by Eqn. [H We also show the eccentricity 
values for their Color Glass Condensate calculation for 
comparison, which interestingly track our x — 1.00 case, 

1. e. follow the density of the binary collisions, except in 
the most central collisions. 

We have also checked the eccentricity fluctuations for 
the different x assumptions to see if this might offer an 
experimental method to discriminate between different x 
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FIG. 1: Glauber Monte Carlo eccentricities. Upper panel 
are the values for e par t{2} as a function of the number of 
participants (calculated with different assumptions). Lower 
panel is the ratio of values relative to the Monte Carlo Glauber 
with £=0.0 case. 



values. The results show that the event-to-event fluctua- 
tions in eccentricity are the same for all x cases within 3- 
6%. As one might expect, the fluctuations are somewhat 
smaller for the x = 1.00 case, since there are a larger num- 
ber of binary collisions to smooth out the fluctuations 
compared with participants. However, this difference is 
affected by the spatial correlations between binary colli- 
sions and might even be further modified with inclusion 
of a nucleon-nucleon hard core potential (i.e. a non-zero 
value of d) . 

In addition to the differences in eccentricities even for 
the case of identical x values, there are a number of dif- 
ferences we have uncovered that we describe below. 



• The average overlap area St is calculated using the 
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following equation: 

St = ±*spfrl - °l y (10) 

corresponding to the area of the tilted ellipse. It 
should be noted that, in the literature, this equa- 
tion appears in with various prefactors (e.g. n [25j |. 
2-7T and in The prefactor An gives a 

somewhat more intuitive result, as it yields the 
correct answer for the simple case of a disc with 
uniform density. This convention is utilized in the 
remainder of this paper. However, since the par- 
ton density is not uniform, it is possible that one 
should consider a smaller area (for example in a 
core/corona type picture). 

In the extraction of the product (<rc s ) from Eqn. [2l 
any multiplicative scaling of the vertical-axis quan- 
tity does not change the extracted results, but 
a multiplicative scaling of the horizontal-axis ~ 
3r 4f" resur ts in a linear change in (<jc s ). However, 
we later note that simple multiplicative scaling of 
these quantities cancels in the determination 77/s. 

Finally, Ref. [l| does not calculate the overlap area 
St in the tilted (i.e. fluctuating event-by-event) 
frame. This turns out to be a minor (1-2%) differ- 
ence, since the cross term a xy turns out to be rela- 
tively small compared to the other two moments. 

• The experiments do not measure the parton den- 
sity at early times, but instead measure dN c h/drj 
after freeze-out. The rapidity density dN c h/dy is 
extracted from the measured dN c h/dr) by multiply- 
ing by a simple estimate of the Jacobian factor of 
1.15. However, this factor depends on the particle 
mix and momentum spectra, which are not mea- 
sured by PHOBOS over their full acceptance. In 
the calculations of Drescher et at, they have con- 
verted the dNch/dr] to dN c h/dy with a factor 1.25; 
however, in our calculations we have used the fac- 
tor utilized by PHOBOS of 1.15[25]. One can then 
convert from charged hadrons to all hadrons by as- 
suming a multiplicative factor of 1.5. It should be 
noted that equating this hadron density to the par- 
ton density requires the additional assumption of 
local parton-hadron duality, with an exact factor of 
1.0 scaling. Each of these factors enters as a linear 
scaling of the horizontal-axis quantity ~ i^p^-- 

• Ref. [l[ calculates the average of the product (i.e. 
(1/St x dN/dy)) in the Monte Carlo approachm, 
for both CGC and Glauber. In contrast, since 
PHOBOS uses the measured (dN/dy) and calcu- 
late (St) for each centrality bin, they determine 
the product of averages (i.e. (1/<St) (dN/dy).). 
If one wants to use the average of the product, 
one cannot directly utilize the PHOBOS measured 



dN c h/dr] as a function of N par t, since one no 
longer has a model of the correlated fluctuations 
between the two quantities. Instead one must use 
a parametrization based on the MCG calculation 
event-by-event. Throughout this paper, we cal- 
culate the average of the product. We have also 
checked in the MCG calculation that taking the 
product of the averages and find differences in the 
extracted fit parameters of order 20%. 

In our calculations, even though e and St are cal- 
culated with different x values, x = 0.13 is always 
used to determine the charged particle multiplicity. 
This is because any case considered should have 
the correct particle multiplicity input (which is not 
described by significantly different x values within 
the two-component model). Thus, one should think 
of varying x as an arbitrary way of modifying the 
initial geometry while maintaining the correct mul- 
tiplicity at a given centrality. If one were certain 
that the charged particle multiplicity constrains the 
geometry only around x ~ 0.13, then these other 
x scenarios might be ruled out. It is also notable 
that for the Cu+Cu data, there is a disagreement 
between this parametrization of the charged par- 
ticle multiplicity at mid-rapidity by approximately 
10% in the most central events. 

One final comment about the experimental uncertain- 
ties, and how they enter the fits. The PHOBOS data are 
presented graphically with vertical lines for the quadra- 
ture sum of statistical and systematic uncertainties at the 
90% Confidence Level. Thus, they are scaled down by 1.6 
to convert them into one standard deviation uncertainties 
for input to a x 2 fit- It is also not clear that the system- 
atic uncertainties included here are uncorrelated point- 
to-point, which raises some issues about their inclusion 
in the standard x 2 fit. Were the correlations known, the 
modified \ 2 fitting procedure developed in (30| could be 
used to properly account for the correlations when deter- 
mining the errors on the fitted parameters. We make an 
estimate of this possible systematic uncertainty correla- 
tion in the next section. 

EXTRACTING 77/s 

Fig. shows the results from our calculations follow- 
ing the previously described prescription, using our MCG 
results with a;=0.00, x=0.13, x=1.00. Given the ex- 
perimental data for v-ij 'e par t{2} as a function of trans- 
verse areal density, the fit procedure determines values 
and uncertainties for (ac s ) and (—)ih- Since only the 
product ac s appears in the fit function, determination of 
the cross-section requires an explicit assumption for the 
sound speed. Taking c s = 1/V3 for definiteness one ob- 
tains a cross-section value (in millibarns). The extracted 
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0.3 



0.2 



0.1 



(a x c s ) = (3.11 x 1A/3) ± (0.48 x 1A/3) [mb] 
(v„/e), u = 0.37 ± 0.03 

v 2 'ih 

Fit total x 2 = 5.16 (with 22 data points) 
T|/s x 4ji= 3.99 ± 0.84 (central , 
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(o x c s ) = (3.96 x 1A/3) ± (0.55 x 1A/3) [mb] 
(V.M.. = 0.30 ± 0.02 

1 2 'ih 

Fit total x 2 = 6-09 (with 22 data points) 
T|/s x 47i= 2.59 ± 0.53 (central AuAu) 
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(a x c s ) = (4.52 x 1A/3) ± (0.64 x 1A/3) [mb] 
K/e). u = 0.24 ± 0.01 

v 2 'ih 

Fit total x 2 = 8-79 (with 22 data points) 
T)/s x 4ji= 1 .78 ± 0.38 (central Au Au) 
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FIG. 2: Shown are the PHOBOS data for v 2 /e par t{2} ver- 
sus (g^r^r) fm~ 2 where the e par t{2} and St are calculated 
from a Monte Carlo Glauber. The figures are using x=0.00, 
£=0.13, a;=1.00 (in the top, middle, and bottom panel respec- 
tively) . Also shown are the best fit results in each case for the 
two fit parameters (a xc s ) and («2/e)i/i. Additionally, follow- 
ing Eqn. 1131 we calculate the rj/s x 4n for each x assumed 
case. The quoted parameter uncertainties are one standard 
deviation fit uncertainties only, assuming the experimental 
systematic uncertainties are point-to-point uncorrelated. 



values are shown in Fig. [2] 

The extracted cross-section can be converted into a 
shear viscosity by using the results of a classical calcu- 
lation [3l| for massless particles at temperature T inter- 
acting with isotropic cross-section a 



r) = 1.264- 



giving 



- = 1.264— 

s sa 



1.264- 



T 
Ana 



0.316 



T 



K) dN 
RS T d V 



(11) 



(12) 



In the second step the entropy was assumed to be given 
by s = An, which follows from a definition of the number 
density using an ideal gas equation of state n = P/T and 
the assumption of massless quanta with pressure P = \e 
(here only e is the energy density) . We note that there is 
already a 10% ambiguity here since it is common practice 
in astrophysics to compute the number densities directly 
from the Bose distribution functions fiij], leading to 3.6 
units of entropy per photon. 

Taking T = 200 MeV by assumption (another num- 
ber with substantial uncertainty) and values for R from 
the Monte Carlo Glauber, one can calculate ri/s for the 
most central Au+Au reactions at ^/snn — 200 GeV. The 
resulting values are shown in each panel of Fig. [2l and 
are summarized with later results in Fig. [5] It is notable 
that the initial conditions with larger eccentricity (e.g. 
x = 1.00) yield smaller values of rj/s. This is counter- 
intuitive since a larger initial eccentricity should require 
a larger viscosity to reduce the «2 down to the experi- 
mentally measured value. This result highlights that the 
extracted value of rj/s comes from the curvature of the 
data points as a function of centrality and the x = 1.00 
case indicates the largest degree of flattening (as seen in 
Fig. |g . 

In principle each parameter has its own physics con- 
tent, but it is interesting to check which parameters 
are necessary in order to extract 77/s directly. Since 
K = naR, one can directly write: 



- = 0.32— = 0.32K TR 
s na 



1 



(13) 



and quantities such as the speed of sound c s and the 
transverse overlap area St do not explicitly appear. As 
a result, this form has the considerable advantage of ex- 
plicitly demonstrating the parametric dependence of 77/s 
on the input parameters. Conversely, this simple expres- 
sion holds only for the 'minimalist' assumption for the 
if-dependence given by Eqn. [2l and therefore imposes 
the need to investigate alternative parametrizations of 
this dependence, as discussed in the following section. 

It is interesting to compare our results with those from 
Ref. [H, who obtain rj/s x An = 2.38 using Monte Carlo 
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Glauber x — 0.20 (with the previously noted caveat re- 
garding their non-standard x definition) and 77/s x 4n — 
1.38 using Color Glass Condensate initial conditions. Be- 
fore comparing these values to our results, we note that 
there is a mistake in their calculation of 77/s in the CGC 
case. In the CGC case, in extracting the cross-section a, 
they assume a speed of sound c s = 0.73 x l/v3) but then 
use a density n calculated with the input c s = 1/V3- 
As can be seen explicitly in our Eqn. [T3l the speed of 
sound cancels out. Thus, the correct CGC value from 
their analysis is 77/s x Att — 1.89. We also note that their 
quoted value of n = 3.9 fm~ 3 from reference [3] does 
not appear in that Letter, and must be converted for the 
difference in the St prefactor, the use of optical Glauber 
versus Monte Carlo differences, and a missing Jacobian 
factor. 

Taking these observations into account, we note that 
the values from our calculation for x = 0.13 with 77/s x 
4n = 2.59±0.53 for central Au+Au is similar to that from 
Ref. [H for their x = 0.20 with 77/s x 4vr = 2.38 (despite 
the mismatch in exact eccentricities and other differences 
outlined above). Similarly, our calculation with x = 1.00 
gives 77/s x 4ir — 1.78 ± 0.38, similar to Ref. [l| for the 
CGC case with 77/s x An — 1.89 (with the correction 
described above). 

Here we re- visit the issue of the handling of the PHO- 
BOS experimental systematic uncertainties. In the above 
we have considered them to be point-to-point uncorre- 
cted. Since we do not know the full correlation matrix, 
we consider two cases to estimate of the systematic uncer- 
tainty on the extracted 77/s. If we allow the data points 
to move within the one standard deviation uncertain- 
ties in a correlated or anti-correlated manner (i.e. tilt- 
ing the i>2/epari{2} values around the mid-central point), 
we find that for the x = 0.13 case, we obtain values of 
77/s = 1.68 (moving the central points down and the pe- 
ripheral points up) and 77/s — 3.61 (moving the central 
points up and the peripheral points down) . We note that 
the x 2 total is very large in all cases in part because the 
statistical uncertainties are quite small and some part of 
the systematic uncertainty is likely uncorrelated. Thus, 
as an approximate estimate, the true fit uncertainty from 
the experimental uncertainties alone is most likely twice 
as large as that shown in Fig. [2j For the x — 0.13 case, 
we then should state it as 77/s = 2.59 ± 1.00, instead of 
the previously quoted 77/s — 2.59 ± 0.53. 



OTHER PARAMETRIZATIONS 

It has already been noted that Eqn. [2] is not derived 
from any a priori theoretical expectation, but has been 
constructed to obey limiting functional forms at large 
and small K. Given this, it is of interest whether an 
alternative functional form might fit the available data 
equally well and still obey the same two limits. A sim- 
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(a x c s ) = (4.97 x 1A/3) ± (2.38 x 1A/3) [mb] 
(v„/e), u = 0.22 ± 0.02 

v 2 'ih 

Fit total x 2 = 5.46 (with 22 data points) 
t|/s x 4ji= 2.04 ± 0.79 (central AuAu) 
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FIG. 3: Shown are the PHOBOS data for V2/e par t{'2} versus 
(s^^r) f m 2 where the eccentricity and transverse overlap 
area are calculated from a Monte Carlo Glauber. The calcu- 
lation uses x=0.13. Also shown are the best fit results (solid 
line) using the formulation with the Pade approximant moti- 
vated functional form, and then extracting the fit parameters 
(a x c 3 ) and vijtih). Additionally, following Eqn. 1121 we cal- 
culate the 77/s x 4-7T. For comparison, we plot the curve for 
the same case with the standard fit (dashed line). 



pie functional form (motivated by Pade approximants) is 
utilized here: 



l + A 



K 

A'n 



D 



/// 



1 



(14) 



This equation obeys the large and small K limits as 
the previous parametrization, if B = E and C = A + 1. 
Of course one can obtain an infinite set of such equations 
by expanding the number of terms in the numerator and 
denominator. For this case, we constrain ourselves to 
the case where B = E = 1 and A = 1, C = 2. Fig. [3] 
shows the best fit where the additional parameter D is 
allowed to vary (for the x = 0.13 case). There is an 
approximately 25% decrease in the ideal hydrodynamic 
value (^f-)ih, while the a = 4.97 ± 2.38 millibarns (as- 
suming c s = l/v3) increases by about 25% relative to 
the previous value of 3.96 ± 0.55, although with a larger 
uncertainty. This best fit corresponds to D sa 2. 

However, it turns out that the standard MINUIT \ 2 
fit has returned an incorrect uncertainty range because 
of a non-parabolic shape of the \ 2 surface. The reason 
for this behavior is examined in Fig. [H which shows all 
of the parameter fits with \ 2 values within one standard 
deviation of the minimum, the total \ 2 as a function of 
the extracted o~c s value, and the associated parameter 
values for D and — ., . 

e %h 
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to determine the functional form, then one needs to make 
a proper theoretical argument about which types func- 
tional forms are truly relevant. There have been attempts 
to check the Knudsen number dependence with trans por t 
calculations and viscous hydrodynamic calculations [33j |. 
and these may lend support to the simple parametriza- 
tion. However, if the parametrization is only tested in 
certain scenarios, then one has all the limitations of those 
particular scenarios and one cannot make more general 
conclusions. 



QUANTUM LIMITS 



(cxc,)/(1Al3)[mb] 



(<rxc,)/(1^)[mb] 



FIG. 4: Upper panel shows the fit results with \ 2 values 
within la of the minimum. Middle panel shows the total \ 2 
as a function of the extracted ac 3 value. Lower panel shows 
the associated parameter values for D (dashed) and 2* xlO 
for visibility (solid). 



The presented % 2 as a function of (ac s ) allows the de- 
termination of the uncertainty on a at one-standard de- 
viation to be ^^ 5 q° millibarns. Almost any arbitrarily 
large value for the cross-section a gives an equally good 
description of the experimental data with an almost iden- 
tical value for — . , and an increasing value of D for larger 
a values. Using Eqn. [12] to determine the corresponding 
values for r)/s we find an allowed range of rj/s x A-k from 
2.55 all the way down 0.34 (well below the postulated 
bound) . 

These results indicate a certain fragility in the proce- 
dure outlined in Ref. [1[ since a minimal modification of 
the parametrization amplifies the uncertainties in rj/s by 
an order of magnitude. It is therefore important to inves- 
tigate if there are a priori theoretical restrictions on the 
functional dependence of v^jt on the Knudsen number. 
For example, are there theoretical arguments that might 
constrain D ~ 2. In a simple transport scenario with cen- 
trality independent cross-section a, temperature T, and 
speed of sound c s , one might question why the underly- 
ing dynamics would change significantly with centrality 
(i.e. D » 1). 

However, if obeying the two limits for K is insufficient 



One obvious feature of the above treatment is the lack 
of any explicit quantum bound on the value for i]/s. 
Clearly if such a bound exists, then the experimental 
values for v-ijt cannot approach the ideal hydrodynamic 
limit, and it is unlikely that a transport formalism that 
allows violations of the bound would be applicable to 
the determination of rj/s in the vicinity of the bound. In 
this section we consider a modification of the formulation 
which explicitly incorporates the bound. 

Following the original formulation of the bound [3], 
we assume that the mean free path A cannot be smaller 
than the DeBroglie wavelength of the particle. A sim- 
ple prescription incorporating this limit is to modify the 
mean free path accordingly: 



no no 



1 



1 



; -<p>A 



(15) 



In the limit where (p) <C n<r, A = l/(p) (the DeBroglie 
wavelength) , while the limit (p) ^> no gives the standard 
expression A = l/(na). Note that we are using natural 
units throughout. 

In relating these values to rj/s we proceed as before 



^ = 0.316^ 

s no 



0.316 



T 



na\l - e-^/ na \ 



(16) 



Taking (p) w 2.7T (relevant for massless, non-interacting 
particles) the high density limit gives rj/s w 0.316 Tx A = 
0.316T/(2.7T) = Thus, there is a modest inconsis- 
tency (at the 40% level) between this implementation and 
the exact bound value. Note that it is possible to enforce 
the precise value of the bound by requiring A > 0.7/(p), 
but in the interest of a simpler heuristic treatment we do 
not consider this additional modification. We highlight 
that even in the presence of quantum effects the Knud- 
sen number itself K = A/ R can come arbitrarily close to 
zero, for example in a neutron star where R ^> any other 
scale in the problem. 



The expression for A' as a function of 



1 dN 

St dy 



is: 



R Rna 



1 



I _ e -2.7T/ncr 



(17) 
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(a x c s ) = (3.45 x 1A/3) ± (0.43 x 1A/3) [mb] 
(v/e) lu = 0.34 ± 0.02 

v 2 'ih 

Fit total x 2 = 6.33 (with 22 data points) 
T|/s x 4ji= 3.43 ±1.17 (central AuAu} 
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FIG. 5: Shown are the PHOBOS data for V2/e par t{'2} versus 
(s^^jf) f m_2 where the eccentricity and transverse overlap 
area are calculated from a Monte Carlo Glauber. The calcula- 
tion uses x=0.13. Also shown are the best fit results using the 
formulation with the quantum bound correction limit (solid 
line), and then extracting the two fit parameters (a x c s ) and 
(v2/e)ih- Additionally, following Eqn. I12I we calculate the 
r\l s x 4-7T. For comparison, we plot the curve for the same 
case without the quantum bound correction (dashed line). 



K 



1 



1 - e 



-2.7T,R/(crc s (; 



dy l> 



(18) 



However, when using this expression into Eqn. [5] one 
needs to know the dependence of R on i^-g^^j before 
performing the fit. This dependence has been calculated 
with a Monte Carlo Glauber and then included in the fit 
to the experimental data. 

The resulting value of the fit is shown in Fig. It gives 
a slightly smaller value of ac s = 3.45±0.43(xl/-\/3) mil- 
libarns. This modest change is not surprising since the 
starting value without the correction is already a factor 
of 2.5 away from the quantum bound. For systems even 
closer to the bound, the impact would be much more sig- 
nificant. Note however that the extracted value of 77/s 
is almost 1 /4-7T larger than that extracted without incor- 
porating the bound. This is not at all unexpected since 
one has re- interpreted the associated 77/s value with the 
asymptotic limit of the fit function. 



DISCUSSION 

In this paper, we are not advocating the Knudsen 
number formalism as a rigorous method of extracting 
77/s. Rather we are highlighting various issues with this 
methodology and note areas where systematic uncertain- 
ties are common to various methods that have been pro- 
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FIG. 6: The solid points are the extracted values of 77/s as 
a function of Monte Carlo Glauber parameter x and for the 
quantum limit and modified Knudsen number parametriza- 
tion cases. The vertical uncertainty lines are the one standard 
deviation uncertainties from the fits to the experimental data 
assuming the PHOBOS systematic uncertainties are point-to- 
point uncorrelated. An estimate of the systematic uncertainty 
on ?7/s, if the experimental systematic uncertainties are cor- 
related (as discussed in the text), is shown as a box for the 
x = 0.13 case. Also shown as open symbols are the values 
extracted by Drescher et al. for Glauber and CGC initial 
conditions. 



posed for extracting 77/s. There are three main conclu- 
sions. 

1. As has been previously observed, any method for 
extracting 77/s will be sensitive to the details of the initial 
conditions, and in particular the spatial distribution (and 
fluctuations) of the deposited energy density. This in- 
duces a strong sensitivity to the assumed ratio x of binary 
collisions to participants in the Monte Carlo Glauber for- 
malism since a modest ccntrality dependent change in 
the spatial distribution substantially changes the asymp- 
totic behavior implied by the fit. This is particularly 
true because the experimental data does not have small 
enough uncertainties for central Au+Au events to deter- 
mine whether an asymptotic limit has been reached. This 
is the main reason why choosing a different functional 
form for the Knudsen number parametrization allows al- 
most any value for 77/s to be extracted, as demonstrated 
above. 

In this paper we do not consider the Color Glass Con- 
densate initial conditions in detail. However we note 
that there are at least two formulations that give signif- 
icantly different results. These include calculations such 
as Ref. [35| which do not include fluctuations, and Ref. 
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[261 ] which includes fluctuations and is the basis of the 
calculations by Drescher et al. shown in Fig. [T] The 
differences in CGC initial condition eccentricities are as 
large as the differences between Monte Carlo Glauber x 
values, and thus this also needs to be addressed and rec- 
onciled. 

These initial geometry uncertainties generally have a 
larger influence on this extraction method since a subtle 
difference in centrality dependence (i.e. curvature) can 
dramatically change the fit. Note that while Ref. [25I ] 
finds the change in eccentricity from x — 0.0 to x = 0.13 
for central Au+Au events to be 15% and thus considered 
small, we find that the change in the centrality depen- 
dence modifies the extracted value of 77/s by more than 
50%. 

These uncertainties certainly also influence other 
methods of extracting r\j s. We note that in the viscous 
hydrodynamic calculations in Ref. [12J , they consider two 
sets of initial conditions. In the Glauber case they uti- 
lize an optical model without fluctuations, and utilize 
the x = 1.00 case. In contrast, Ref. [10| utilizes partici- 
pants for initial conditions (i.e. x = 0.00). In the case of 
constraints on 77/s from charm and bottom elliptic flow 
V2 @ the extractions of 77/s will also depend on the initial 
conditions and the same scrutiny and propagation of un- 
certainties is needed. In Ref. [36J they utilize an optical 
Glauber model without fluctuations and wounded nuclc- 
ons (i.e. x — 0.00), while in Ref. [37| they state that for 
impact parameter b = 7 fm the eccentricity is e = 0.6. 

These observations highlight the critical importance of 
a reproducible set of publicly available code for generat- 
ing what are often simply referred to as "Glauber" initial 
conditions. All publications need a detailed specification 
of all relevant features of the Glauber calculation, in- 
cluding 1) the x value, 2) Glauber parameters for the 
Woods-Saxon distribution, 3) any d exclusion radius and 
re- weighting of the radius value, 4) any outer radius cut- 
off, 5) whether fluctuations in the reaction plane direction 
are included (as seen in Eqn. [9]) 6) whether the result is 
not a Monte Carlo but instead an optical Glauber ex- 
traction. Many of the individual uncertainties related to 
these parameters have previously been considered minor 
(i.e. 20% or less). 

2. Even if one accepts the assumptions of the formal- 
ism from Ref. [l|, there are significant uncertainties in 
the parameters characterizing the medium on time scales 
corresponding to the largest gradients, i.e. T, c s and 
R. Since these values vary as a function of both the 
space and time coordinates, one needs to specify the pre- 
cise time and space interval over which the averaging is 
performed. It is also necessary to demonstate that such 
averaging of intermediate quantities is handled in a self- 
consistent fashion when deriving other averaged quanti- 
ties (such as r)/s). For example, the speed of sound c s is 
introduced into Eqn. [5] when estimating the characteris- 
tic time over which gradients develop. Absent a detailed 



simulation, the speed of sound in such an expression is 
an ill-determined average over times prior to the desired 
characteristic time. Nor is it obvious that this is same 
time-average assumed for other quantities such as T, R, 
A, etc. 

Moreover, there is also a potential inconsistency with 
the parameters used in studies up to this point. The tem- 
perature, parton density, speed of sound, and interaction 
cross-sections are not independent parameters. In fact 
at zero baryon density, the temperature determines all 
of the others. As an explicit example, a temperature of 
T = 200 MeV is inconsistent with the speed of sound of 
1/a/3 when comparing with lattice QCD results at zero 
baryon chemical potential [38j |. Also, temperature and 
parton density can be translated into T and e (energy 
density) via the QCD equation of state. While it seems 
reasonable that the temperature, density and speed of 
sound are approximately constant over the range of den- 
sities for Cu+Cu and Au+Au collisions at a single col- 
liding energy (e.g. ^/snn = 200 GeV) considered here, 
this is not expected to be true at asymptotically high 
densities. Even though there are no experimental data 
points at such densities, and one extracts a x c s and 77/s 
at the density corresponding to Au+Au central collisions, 
even small changes can systematically alter the central- 
ity dependence, modifying the extracted value of K, and 
thus 77/s. One must also consider the possibility of tem- 
perature dependent cross-sections in medium, something 
not addressed so far. Again, what are often dismissed as 
small centrality dependencies can modify the curvature 
of the data and have an unexpectedly large impact on 
77/s. 

We note that there has been a recent attempt to cali- 
brate the simple analytic approach discussed here, utiliz- 
ing viscous hydrodynamics calculations [33| . These stud- 
ies have already revealed an ambiguity in the value of 
77/s and the equation of state. Additionally, any such 
calibration will have the same caveats and limitations of 
the existing viscous hydrodynamic calculations (with one 
such example being the absence of a hadronic phase after 
freeze-out). 

3. In the formulation discussed in this paper, the ex- 
tracted value of 77/s depends on whether or not a lower 
bound is explicitly introduced. Thus, if a lower bound 
on the viscosity does in fact exist, then 77/s extractions 
which do not incorporate it in their formalism will ulti- 
mately be unreliable. Conversely, approaches which en- 
force such a quantum bound on 77/s obviously cannot ex- 
tract values below the bound. Thus one could never fal- 
sify the existence of such a bound. This is in contrast to 
full viscous hydrodynamic calculations and comparisons 
which do not rely on microscopic modeling of transport 
phenomena. 

This raises the question at to whether the calculational 
framework discussed in Drescher et al, and this work, is 
valid in the strong-coupling limit, i.e. at 77/s = 1/A-k. A 
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standard Boltzmann transport calculation is the oppo- 
site limit from hydrodynamic calculation, which assume 
a continuum hypothesis. The theory community has ex- 
plored both approaches with great benefit, and in both 
cases one can attempt to push the models past their nom- 
inal domains of applicability. If the viscous corrections 
to the ideal hydrodynamics are too large, then they vio- 
late the gradient expansion made in the small rj/s limit. 
In the Boltzmann transport approach (e.g. a parton cas- 
cade), the quasi-particle widths become of same order as 
the masses when the mean free paths are comparable to 
the DeBroglie wavelengths (i.e. T/m R- 1) |___|. As an 
example, if one ignored multi-particle quantum effects in 
a gas of cold atoms near the BEC transition that would 
have a dramatic impact on the physics picture. Thus, in 
this case, it is not clear how to determine the systematic 
uncertainties and overall physical picture. 



SUMMARY 

The formalism for relating rj/s to the Knudsen num- 
ber has attracted much recent attention since it offers a 
tantalizing possibility of extracting transport properties 
almost directly from experimental data. However, several 
important issues need to be addressed consistently before 
one can consider this approach to produce anything more 
than order-of-magnitude estimates. The linearity of rj/s 
with T, c s , and even R makes a serious evaluation of 
their uncertainties of utmost important. One must also 
consider the fact that c s is a function of temperature due 
to the QCD equation of state itself. The final result for 
rj/s is quite sensitive to Monte Carlo Glauber parame- 
ters, and particularly the parameter x which controls the 
balance between N part and N co u. Finally, implementing 
a lower bound to rj/s, for instance as we have presented 
in this work, has a substantial effect on the extracted 
value, particularly since previous calculations find values 
so close to the bound already. 
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